﻿#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cuda.h"
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <Windows.h>
#include <time.h>

#define HANDLE_ERROR(err) if (err != cudaSuccess){\
	fprintf(stderr, "%s: [%s , %d]\n", cudaGetErrorString(err), \
	__FILE__, __LINE__); exit(EXIT_FAILURE); }

#define HANDLE_NULL(ptr) if (!ptr){\
	fprintf(stderr, #ptr" = NULL: [%s , %d]\n", __FILE__, \
	__LINE__); \
	exit(EXIT_FAILURE); }

#define INDEX(x,y,stride,channels) ((y)*(stride) + (x)*(channels))
#define CUBE(x) (x*x*x)
#define SQUARE(x) (x*x)
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))

typedef unsigned char uchar;
typedef unsigned int uint;

typedef struct tagPoint2f
{
	double x;
	double y;
	tagPoint2f(){
		x = 0, y = 0;
	}
	tagPoint2f(double _x, double _y){
		x = _x, y = _y;
	}
} Point2f;
typedef struct tagBlockData
{
	int idu, idv, rows, cols, channels;
	BYTE* dataPtr;
	tagBlockData(){
		idu = idv = rows = cols = channels = 0;
		dataPtr = NULL;
	}
	tagBlockData(int _idu, int _idv, int _rows, int _cols,
		int _channels, BYTE *_dataPtr){
		idu = _idu, idv = _idv;
		rows = _rows, cols = _cols;
		channels = _channels;
		dataPtr = _dataPtr;
	}
} BlockData;
std::vector<BlockData> queueData;

enum
{
	NEARESTNEIGHBOR,
	BILINEAR,
	BICUBIC,
};

__device__ const double X_1_dev[4][4]
{
	{ -1.0 / 6, 1.0 / 2, -1.0 / 3, 0.0 },
	{ 1.0 / 2, -1.0, -1.0 / 2, 1.0 },
	{ -1.0 / 2, 1.0 / 2, 1.0, 0.0 },
	{ 1.0 / 6, 0.0, -1.0 / 6, 0.0 }
};

__device__ const double Y_1_dev[4][4]
{
	{ -1.0 / 6, 1.0 / 3, -1.0 / 2, 1.0 / 6 },
	{ 1.0 / 2, -1.0 / 2, 1.0 / 2, 0.0 },
	{ -1.0 / 3, -5.0 / 6, 1.0, -1.0 / 6 },
	{ 0.0, 1.0, 0.0, 0.0 }
};

const double X_1_hst[4][4]
{
	{ -1.0 / 6, 1.0 / 2, -1.0 / 3, 0.0 },
	{ 1.0 / 2, -1.0, -1.0 / 2, 1.0 },
	{ -1.0 / 2, 1.0 / 2, 1.0, 0.0 },
	{ 1.0 / 6, 0.0, -1.0 / 6, 0.0 }
};

const double Y_1_hst[4][4]
{
	{ -1.0 / 6, 1.0 / 3, -1.0 / 2, 1.0 / 6 },
	{ 1.0 / 2, -1.0 / 2, 1.0 / 2, 0.0 },
	{ -1.0 / 3, -5.0 / 6, 1.0, -1.0 / 6 },
	{ 0.0, 1.0, 0.0, 0.0 }
};

__device__ void NearestNeighbor_dev(const uchar *inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v);

__device__ void Bilinear_dev(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v);

__device__ void Bicubic_dev(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v);

void NearestNeighbor_hst(const uchar *inPtr, int strideIn,
	int channels, int x, int y,
	uchar *outPtr, int strideOut, double u, double v);

void Bilinear_hst(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v);

void Bicubic_hst(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v);

__global__ void TransformKernel(
	uchar *outPtr, int u0, int v0, 
	int outHeight, int outWidth, int channels,
	const double *g,
	const uchar *inPtr, int x0, int y0, 
	int inHeight, int inWidth,
	int mode = 0);

void DocMauAnh(FILE* filePtr, BYTE headerSize,
	int r0, int c0, int rows, int cols, int bytesPerPixel,
	int rowLength, BYTE* &outPtr);

bool DocPhanHeader(FILE *filePtr, int &width, int &height, int &depth,
	int &rowLength, int &headerSize);

void Inv(const double f[6], double g[6]);

void TransformWithCPU(BYTE *&outPtr, int u0, int v0,
	int outWidth, int outHeight, int channels, const double *g,
	const BYTE *inPtr, int x0, int y0,
	int inWidth, int inHeight, int mode);

void TransformWithGPU(BYTE *&outPtr, int u0, int v0,
	int outWidth, int outHeight, int channels, const double *g,
	const BYTE *inPtr, int x0, int y0,
	int inWidth, int inHeight, int mode,
	dim3 threadsPerBlock = dim3(4, 4),
	dim3 blocksPerGrid = dim3(4, 4));

void GetBox(int u0, int v0, int rows, int cols, const double *f,
	double &minx, double &miny, double &maxx, double &maxy);

void help();



int main(int argc, char* argv[])
{
	int i = 0, k = 0;
	int iRows = 512; //Số hàng pixel của một mẩu ảnh kết quả
	int iCols = 512; //Số pixel trên một hàng của một mẩu ảnh kết quả
	double f[6] = { 1, 0, 0, 0, 1, 0 }; // hàm biến đổi
	char processor[6] = "cpu"; // chế độ chạy CPU | GPU | OpenMP
	int mode = 0; // phương pháp tái chế mẫu NEARESTNEIGHBOR | BILINEAR | BICUBIC
	dim3 threadsPerBlocks(4, 4); // số threads 2D trên một block
	dim3 blocksPerGrid(4, 4); // số blocks 2D trên grid

	//---- Lấy các tham số biến đổi ----
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/f") == 0){
			if (i + 6 < argc)
				for (k = i + 1; k < i + 7; ++k)
					f[k - i - 1] = atof(argv[k]);
			else{
				help();
				return 1;
			}
			break;
		}
	}

	//---- Lấy phương pháp tái chế mẫu ----
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/m") == 0){
			if (i + 1 < argc){
				if (strcmp(strlwr(argv[i + 1]), "bilinear") == 0)
					mode = 1;
				else if (strcmp(strlwr(argv[i + 1]), "bicubic") == 0)
					mode = 2;
				else mode = 0;
			}
			else{
				help();
				return 1;
			}
			break;
		}
	}

	//---- Lấy chế độ sẽ chạy ----
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/p") == 0){
			if (i + 1 < argc){
				if (strcmp(strlwr(argv[i + 1]), "openmp") == 0)
					strcpy(processor, "openmp");
				else if (strcmp(strlwr(argv[i + 1]), "gpu") == 0)
					strcpy(processor, "gpu");
				else strcpy(processor, "cpu");
			}
			else{
				help();
				return 1;
			}
			break;
		}
	}

	//---- Mẩu ảnh sẽ trích để chạy trong một lần ----
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/s") == 0){
			if (i + 2 < argc){
				iRows = atoi(argv[i + 1]);
				iCols = atoi(argv[i + 2]);
			}
			else{
				help();
				return 1;
			}
			break;
		}
	}

	//---- Số lượng threads 2D trên một blocks ----
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/t") == 0){
			if (i + 2 < argc){
				threadsPerBlocks = dim3(
					atoi(argv[i + 1]), atoi(argv[i + 2]));
			}
			else{
				help();
				return 1;
			}
			break;
		}
	}

	//---- Số lượng blocks 2D trên grid
	for (i = 1; i < argc; ++i){
		if (strcmp(strlwr(argv[i]), "/b") == 0){
			if (i + 2 < argc){
				blocksPerGrid = dim3(
					atoi(argv[i + 1]), atoi(argv[i + 2]));
			}
			else{
				help();
				return 1;
			}
			break;
		}
	}



	//---- Đọc ảnh đầu vào ----
	FILE *filePtr = fopen(argv[1], "r"); // con trỏ đọc ảnh đầu vào
	HANDLE_NULL(filePtr);
	int inWidth = 0;	// chiều rộng của ảnh đầu vào
	int inHeight = 0; // chiều cao của ảnh đầu vào
	int depth = 0;  // độ sâu màu (byte)
	int channels = 0; // số kênh ảnh
	int rowLength = 0; // Số bytes trên một hàng ảnh
	int headerSize = 0; // kích thước phần header (bytes)
	bool success = true;
	success = DocPhanHeader(filePtr, inWidth, inHeight, depth, rowLength,
		headerSize);
	channels = depth >> 3;
	if (success == false){
		fprintf(stderr, "KHONG the doc header file");
		return 1;
	}

	//---- Xác định tạo độ và kích thước của ảnh kết quả ----
	Point2f topLeft(0, 0); //Tọa độ góc trên bên trái của ảnh kết quả
	Point2f bottomRight(0, 0); // Tạo độ góc dưới bên phải của ảnh kết quả
	GetBox(0, 0, inWidth, inHeight, f,
		topLeft.x, topLeft.y,
		bottomRight.x, bottomRight.y);
	int outHeight; // Chiều cao của ảnh kết quả
	int outWidth; // Chiều rộng của ảnh kết quả
	outWidth = (int)(ceil(bottomRight.x) - floor(topLeft.x));
	outHeight = (int)(ceil(bottomRight.y) - floor(topLeft.y));

	//---- Thực hiện biến đổi trên từng mẩu ảnh ----
	BYTE *inPtr = NULL; //Con trỏ đến vùng dữ liệu ảnh sẽ lấy ra để xử lý
	BYTE *outPtr = NULL; //Con trỏ vùng dữ liệu ảnh sau khi đã xử lý
	int rows = 0; // Số hàng pixel sẽ đọc
	int cols = 0; // Số pixel sẽ đọc trên một hàng
	double g[6]; // Tham số của hàm ánh xạ ngược
	Inv(f, g);
	Point2f subTopLeft(0, 0); // Tạo độ góc trên bên trái của mẩu ảnh gốc
	Point2f subBottomRight(0, 0); // Tạo độ góc dưới bên phải của mẩu ảnh gốc
	int subRows = 0; // Số hàng pixel của mẩu ảnh gốc
	int subCols = 0; // Số pixel trên một hàng cua mẩu ảnh gốc
	int all = int(ceil(outHeight*1.0 / iRows)*ceil(outWidth*1.0 / iCols));
	int count = 0;
	outPtr = (BYTE*)malloc((iRows+1)*(iCols+1)*channels*sizeof(BYTE));
	HANDLE_NULL(outPtr);

	if (strcmp(processor, "gpu") == 0){
		cudaEvent_t start, stop;
		float kernelTime = 0;
		for (int v0 = 0; v0 < outHeight; v0 += iRows){
			for (int u0 = 0; u0 < outWidth; u0 += iCols){
				rows = MIN(iRows, outHeight - v0);
				cols = MIN(iCols, outWidth - u0);
				GetBox(u0, v0, cols, rows, g,
					subTopLeft.x, subTopLeft.y,
					subBottomRight.x, subBottomRight.y);
				subCols = (int)(ceil(subBottomRight.x) - floor(subTopLeft.x));
				subRows = (int)(ceil(subBottomRight.y) - floor(subTopLeft.y));
				if (!inPtr){
					inPtr = (BYTE*)malloc((subRows+1)*(subCols+1)*channels*sizeof(BYTE));
					HANDLE_NULL(inPtr);
				}
				//printf("\n%d %d \t %d %d\n", subRows, subCols, rows, cols);
				
				DocMauAnh(filePtr, headerSize,
					(int)subTopLeft.y, (int)subTopLeft.x, rows, cols,
					channels, rowLength	, inPtr);
				cudaEventCreate(&start);
				cudaEventRecord(start, 0);
				TransformWithGPU(outPtr, u0, v0, cols, rows, channels, g,
					inPtr, (int)subTopLeft.x, (int)subTopLeft.y,
					subCols, subRows, mode, threadsPerBlocks, blocksPerGrid);
				cudaEventCreate(&stop);
				cudaEventRecord(stop, 0);
				cudaEventSynchronize(stop);
				cudaEventElapsedTime(&kernelTime, start, stop);
				queueData.push_back(BlockData(u0, v0, rows, cols, channels, outPtr));
				//putchar(13);
				printf("\nDang thuc hien: %d/%d = %0.2f%%", ++count, all, count*100.0 / all);
			}
		}
		printf("\nThoi gian: %0.2fms", kernelTime);
	}
	else if (strcmp(processor, "openmp") == 0){
		printf("Chua cai dat\n");
	}
	else{
		int beginTime = clock();
		for (int v0 = 0; v0 < outHeight; v0 += iRows){
			for (int u0 = 0; u0 < outWidth; u0 += iCols){
				rows = MIN(iRows, outHeight - v0);
				cols = MIN(iCols, outWidth - u0);
				GetBox(u0, v0, cols, rows, g,
					subTopLeft.x, subTopLeft.y,
					subBottomRight.x, subBottomRight.y);
				subCols = (int)(ceil(subBottomRight.x) - floor(subTopLeft.x));
				subRows = (int)(ceil(subBottomRight.y) - floor(subTopLeft.y));
				if (!inPtr){
					inPtr = (BYTE*)malloc(subRows*subCols*channels*sizeof(BYTE));
					HANDLE_NULL(inPtr);
				}
				//printf("\n%d %d \t %d %d\n", subRows, subCols, rows, cols);

				DocMauAnh(filePtr, headerSize,
					(int)subTopLeft.y, (int)subTopLeft.x, rows, cols,
					channels, rowLength, inPtr);

				TransformWithCPU(outPtr, u0, v0, cols, rows, channels, g,
					inPtr, (int)subTopLeft.x, (int)subTopLeft.y,
					subCols, subRows, mode);
				queueData.push_back(BlockData(u0, v0, rows, cols, channels, outPtr));
				putchar(13);
				printf("Dang thuc hien: %d/%d = %0.2f%%", ++count, all, count*100.0 / all);
			}
		}
		printf("\nThoi gian: %0.2fs", float(clock() - beginTime) / CLOCKS_PER_SEC);
	}
	free(inPtr);
	free(outPtr);
	inPtr = outPtr = NULL;
	printf("\nXong !\n");
	fclose(filePtr);
	return 0;	
}



void NearestNeighbor_hst(const uchar *inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	x = int(x + 0.5);
	y = int(y + 0.5);
	int i = INDEX(int(x), int(y), strideIn, channels);
	int j = INDEX(u, v, strideOut, channels);
	for (int r = 0; r < channels; ++r){
		outPtr[j + r] = inPtr[i + r];
	}
}

void Bilinear_hst(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	double dx = x - (int)x;
	double dy = y - (int)y;
	double p1;
	double p2;
	int i = INDEX((int)x, (int)y, strideIn, channels);
	int j = INDEX(u, v, strideOut, channels);
	for (int r = 0; r < channels; ++r){
		p1 = inPtr[i + r] * (1 - dx) + inPtr[i + channels + r] * dx;
		p2 = inPtr[i + strideIn + r] * (1 - dx) +
			inPtr[i + strideIn + channels + r] * dx;
		outPtr[j + r] = uchar(p1*(1 - dy) + p2*dy);
	}
}

void Bicubic_hst(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	double dx = x - (int)x, dy = y - (int)y;
	int i, r, c;
	double temp1[4]{ 0, 0, 0, 0 }; // [Py]*[Y_1_dev]
	for (i = 0; i < 4; ++i){
		temp1[i] = CUBE(dy)*Y_1_hst[0][i] + SQUARE(dy)*Y_1_hst[1][i]
			+ dy*Y_1_hst[2][i] + Y_1_hst[3][i];
	}
	double temp2[4]{0, 0, 0, 0}; // [X_1_dev]*[Px]
	for (i = 0; i < 4; ++i){
		temp2[i] = X_1_hst[i][0] * CUBE(dx) + X_1_hst[i][1] * SQUARE(dx)
			+ X_1_hst[i][2] * dx + X_1_hst[i][3];
	}
	i = INDEX((int)x, (int)y, strideIn, channels);
	uchar* Z[4][4];
	for (r = -1; r < 3; ++r){
		for (c = -1; c < 3; ++c){
			Z[r + 1][c + 1] =
				(uchar*)(inPtr + i + r*strideIn + c*channels);
		}
	}
	double pyz[4][4]; // [Py]*[_Y_1]*Z
	for (r = 0; r < 4; ++r){
		for (c = 0; c < 4; ++c){
			pyz[r][c] = 0;
		}
	}
	for (c = 0; c < channels; ++c){
		for (i = 0; i < 4; ++i){
			pyz[i][c] = temp1[0] * (*(Z[0][i] + c)) +
				temp1[1] * (*(Z[1][i] + c)) +
				temp1[2] * (*(Z[2][i] + c)) +
				temp1[3] * (*(Z[3][i] + c));
		}
	}
	i = INDEX(u, v, strideOut, channels);
	for (c = 0; c < channels; ++c){
		outPtr[i + c] = (uchar)(pyz[0][c] * temp2[0] +
			pyz[1][c] * temp2[1] +
			pyz[2][c] * temp2[2] +
			pyz[3][c] * temp2[3]);
	}
}

__device__ void NearestNeighbor_dev(const uchar *inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	x = int(x + 0.5);
	y = int(y + 0.5);
	int i = INDEX(int(x), int(y), strideIn, channels);
	int j = INDEX(u, v, strideOut, channels);
	for (int r = 0; r < channels; ++r){
		outPtr[j + r] = inPtr[i + r];
	}
}

__device__ void Bilinear_dev(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	double dx = x - (int)x;
	double dy = y - (int)y;
	double p1;
	double p2;
	int i = INDEX((int)x, (int)y, strideIn, channels);
	int j = INDEX(u, v, strideOut, channels);
	for (int r = 0; r < channels; ++r){
		p1 = inPtr[i + r] * (1 - dx) + inPtr[i + channels + r] * dx;
		p2 = inPtr[i + strideIn + r] * (1 - dx) +
			inPtr[i + strideIn + channels + r] * dx;
		outPtr[j + r] = uchar(p1*(1 - dy) + p2*dy);
	}
}

__device__ void Bicubic_dev(const uchar* inPtr, int strideIn,
	int channels, double x, double y,
	uchar *outPtr, int strideOut, int u, int v)
{
	double dx = x - (int)x, dy = y - (int)y;
	int i, r, c;
	double temp1[4]{ 0, 0, 0, 0 }; // [Py]*[Y_1_dev]
	for (i = 0; i < 4; ++i){
		temp1[i] = CUBE(dy)*Y_1_dev[0][i] + SQUARE(dy)*Y_1_dev[1][i]
			+ dy*Y_1_dev[2][i] + Y_1_dev[3][i];
	}
	double temp2[4]{0, 0, 0, 0}; // [X_1_dev]*[Px]
	for (i = 0; i < 4; ++i){
		temp2[i] = X_1_dev[i][0] * CUBE(dx) + X_1_dev[i][1] * SQUARE(dx)
			+ X_1_dev[i][2] * dx + X_1_dev[i][3];
	}
	i = INDEX((int)x, (int)y, strideIn, channels);
	uchar* Z[4][4];
	for (r = -1; r < 3; ++r){
		for (c = -1; c < 3; ++c){
			Z[r + 1][c + 1] =
				(uchar*)(inPtr + i + r*strideIn + c*channels);
		}
	}
	double pyz[4][4]; // [Py]*[_Y_1]*Z
	for (r = 0; r < 4; ++r){
		for (c = 0; c < 4; ++c){
			pyz[r][c] = 0;
		}
	}
	for (c = 0; c < channels; ++c){
		for (i = 0; i < 4; ++i){
			pyz[i][c] = temp1[0] * (*(Z[0][i] + c)) +
				temp1[1] * (*(Z[1][i] + c)) +
				temp1[2] * (*(Z[2][i] + c)) +
				temp1[3] * (*(Z[3][i] + c));
		}
	}
	i = INDEX(u, v, strideOut, channels);
	for (c = 0; c < channels; ++c){
		outPtr[i + c] = (uchar)(pyz[0][c] * temp2[0] +
			pyz[1][c] * temp2[1] +
			pyz[2][c] * temp2[2] +
			pyz[3][c] * temp2[3]);
	}
}


__global__ void TransformKernel(
	uchar *outPtr, int u0, int v0, 
	int outHeight, int outWidth, int channels,
	const double *g,
	const uchar *inPtr, int x0, int y0, 
	int inHeight, int inWidth,
	int mode)
{
	double x;
	double y;
	int u;
	int v;
	if (mode == NEARESTNEIGHBOR){
		v = threadIdx.y + blockIdx.y * blockDim.y;
		while (v < outHeight){
			u = threadIdx.x + blockIdx.x * blockDim.x;
			while (u < outWidth){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					NearestNeighbor_dev(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
				u += gridDim.x * blockDim.x;
			}
			v += gridDim.y * blockDim.y;
		}
	}
	else if (mode == BILINEAR){
		v = threadIdx.y + blockIdx.y * blockDim.y;
		while (v < outHeight){
			u = threadIdx.x + blockIdx.x * blockDim.x;
			while (u < outWidth){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					Bilinear_dev(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
				u += gridDim.x * blockDim.x;
			}
			v += gridDim.y * blockDim.y;
		}
	}
	else{
		v = threadIdx.y + blockIdx.y * blockDim.y;
		while (v < outHeight){
			u = threadIdx.x + blockIdx.x * blockDim.x;
			while (u < outWidth){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					Bicubic_dev(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
				u += gridDim.x * blockDim.x;
			}
			v += gridDim.y * blockDim.y;
		}
	}
}

//------------ Lấy một mẩu ảnh------------
// headerSize: kích thước phần header
// r0: là tọa độ hàng của mẩu
// c0: là tọa độ cột của mẩu
// rows: kích thước hàng
// cols: kích thước cột
// bytesPerPixel: số lượng bytes cho một pixel
// rowLength: số bytes trên một hàng ảnh
// outPtr: con trỏ đầu ra
void DocMauAnh(FILE* filePtr, BYTE headerSize,
	int r0, int c0, int rows, int cols, int bytesPerPixel,
	int rowLength, BYTE* &outPtr){
	HANDLE_NULL(outPtr);
	BYTE* ptr = outPtr;
	//Thực hiện đọc từng hàng ảnh
	fseek(filePtr, headerSize + r0 *rowLength + c0*bytesPerPixel, SEEK_SET);
	for (int i = 0; i < rows; ++i){
		fread(ptr, sizeof(BYTE), cols*bytesPerPixel, filePtr);
		ptr += cols*bytesPerPixel;
		fseek(filePtr, rowLength - cols*bytesPerPixel, SEEK_CUR);
	}
}

bool DocPhanHeader(FILE *filePtr, int &width, int &height, int &depth,
	int &rowLength, int &headerSize){
	//Đọc phần header file
	BITMAPFILEHEADER bmfh;
	int soByteDocDuoc = 0;
	soByteDocDuoc = (int)fread(&bmfh, sizeof(BITMAPFILEHEADER), 1, filePtr);
	if (soByteDocDuoc != 1 ||
		bmfh.bfType != *(WORD*)"BM"){
		fprintf(stderr, "Khong the duoc file, hoac file sai dinh dang.\n");
		fclose(filePtr);
		return false;
	}

	// kích thước phần thông tin ảnh và dữ liệu ảnh
	int infoSize = bmfh.bfOffBits - sizeof(BITMAPFILEHEADER);
	// Đọc phần header ảnh
	BITMAPINFO* pbmi = NULL;
	pbmi = (BITMAPINFO*)malloc(infoSize);
	HANDLE_NULL(pbmi);
	soByteDocDuoc = (int)fread(pbmi, sizeof(BITMAPINFO), 1, filePtr);
	if (soByteDocDuoc != 1){
		fprintf(stderr, "Khong the duoc file, hoac file sai dinh dang.\n");
		fclose(filePtr);
		return false;
	}
	// Tính chiều dài , chiều rộng
	// số lượng bytes trên một hàng ảnh
	width = pbmi->bmiHeader.biWidth;
	height = pbmi->bmiHeader.biHeight;
	depth = pbmi->bmiHeader.biBitCount;
	rowLength = ((width * depth + 31) & ~31) >> 3;
	headerSize = sizeof(BITMAPFILEHEADER)+infoSize;
	return true;
}

void Inv(const double f[6], double g[6])
{
	g[0] = f[4];
	g[1] = -f[1];
	g[2] = f[1] * f[5] - f[2] * f[4];
	g[3] = -f[3];
	g[4] = f[0];
	g[5] = f[2] * f[3] - f[0] * f[5];
	for (int i = 0; i < 6; i++)
		g[i] /= (f[0] * f[4] - f[1] * f[3]);
}

// ---------- Xác định Box của ảnh ----------
// u0: tạo độ cột
// v0: tạo độ hàng của ảnh
// width: chiều rộng, height là chiều cao
// f: là con trỏ các tham số biến đổi
// Hàm trả về giá trị nhỏ nhất và lớn nhất của x và y;
void GetBox(int u0, int v0, int width, int height, const double *f,
	double &minx, double &miny, double &maxx, double &maxy)
{
	double g[6];
	Inv(f, g);
	Point2f goc[4];
	goc[0] = Point2f(
		f[0] * u0 + f[1] * v0 + f[2],
		f[3] * u0 + f[4] * v0 + f[5]);
	goc[1] = Point2f(
		f[0] * (u0 + width) + f[1] * v0 + f[2],
		f[3] * (u0 + width) + f[4] * v0 + f[5]);
	goc[2] = Point2f(
		f[0] * (u0 + width) + f[1] * (v0 + height) + f[2],
		f[3] * (u0 + width) + f[4] * (v0 + height) + f[5]);
	goc[3] = Point2f(
		f[0] * u0 + f[1] * (v0 + height) + f[2],
		f[3] * u0 + f[4] * (v0 + height) + f[5]);
	minx = maxx = goc[0].x;
	miny = maxy = goc[0].y;
	for (int i = 0; i < 4; ++i){
		if (goc[i].x > maxx)
			maxx = goc[i].x;
		if (goc[i].x < minx)
			minx = goc[i].x;
		if (goc[i].y > maxy)
			maxy = goc[i].y;
		if (goc[i].y < miny)
			miny = goc[i].y;
	}
}

//---------- Hàm biến đổi sử dụng GPU ----------
// outPtr: con trỏ của mẩu ảnh kết quả
// u0: Tọa độ theo cột của mẩu ảnh kết quả
// v0: Tạo độ theo hàng của mẩu ảnh kết quả
// outWidth: chiều rộng của mẩu ảnh kết quả
// outHeight: chiều cao của mẩu ảnh kết quả
// g: tham số của hàm ánh xạ ngược
// inPtr: con trỏ của mẩu ảnh gốc
// x0: Tạo độ theo cột của mẩu ảnh gốc
// y0: Tạo độ theo hàng của mẩu ảnh gốc
// inWidth: chiều rộng của mẩu ảnh gốc
// inHeight: chiều cao của mẩu ảnh gốc
// mode: Phương pháp tái chế mẫu
// threadsPerBlock: Số theads 2D trên một Blocks, mặc định (4,4)
// blocksPerGrid: Số blocks 2D trên Grid, mặc định (4,4)
void TransformWithGPU(BYTE *&outPtr, int u0, int v0,
	int outWidth, int outHeight, int channels, const double *g,
	const BYTE *inPtr, int x0, int y0,
	int inWidth, int inHeight, int mode,
	dim3 threadsPerBlock,
	dim3 blocksPerGrid)
{
	HANDLE_NULL(inPtr);
	HANDLE_NULL(outPtr);
	uchar *inPtrDev = NULL; // con trỏ mẩu ảnh gốc trên device
	uchar *outPtrDev = NULL; // con trỏ mẩu ảnh kết quả trên device
	double *gPtrDev = NULL; // tham số của hàm ánh xạ ngược trên device
	// Cấp phát bộ nhớ
	HANDLE_ERROR(cudaMalloc((void**)&inPtrDev,
		inWidth*inHeight*channels));
	HANDLE_ERROR(cudaMalloc((void**)&outPtrDev,
		outWidth*outHeight*channels));
	HANDLE_ERROR(cudaMalloc((void**)&gPtrDev,
		6 * sizeof(double)));
	//Copy dữ liệu từ host sang device
	printf("\n %X %X %d %d %d",
		inPtrDev, inPtr, inWidth, inHeight, channels);
	HANDLE_ERROR(cudaMemcpy(inPtrDev, inPtr,
		inWidth*inHeight*channels,
		cudaMemcpyHostToDevice));
	HANDLE_ERROR(cudaMemcpy(gPtrDev, g,
		6 * sizeof(double),
		cudaMemcpyHostToDevice));
	// chạy nhân
	TransformKernel <<< blocksPerGrid, threadsPerBlock >>> (
		outPtrDev, u0, v0, outHeight, outWidth,
		channels, gPtrDev,
		inPtrDev, x0, y0, inHeight, inWidth,
		mode);

	HANDLE_ERROR(cudaGetLastError());
	HANDLE_ERROR(cudaMemcpy(outPtr, outPtrDev,
		outWidth*outHeight*channels, cudaMemcpyDeviceToHost));
	cudaFree(inPtrDev);
	cudaFree(outPtrDev);
	cudaFree(gPtrDev);
}


void TransformWithCPU(BYTE *&outPtr, int u0, int v0,
	int outWidth, int outHeight, int channels, const double *g,
	const BYTE *inPtr, int x0, int y0,
	int inWidth, int inHeight, int mode)
{
	double x, y;
	if (mode == NEARESTNEIGHBOR){
		for (int v = 0; v < outHeight; ++v){
			for (int u = 0; u < outWidth; ++u){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					NearestNeighbor_hst(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
			}
		}
	}
	else if (mode == BILINEAR){
		for (int v = 0; v < outHeight; ++v){
			for (int u = 0; u < outWidth; ++u){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					Bilinear_hst(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
			}
		}
	}
	else{
		for (int v = 0; v < outHeight; ++v){
			for (int u = 0; u < outWidth; ++u){
				x = g[0] * (u0 + u) + g[1] * (v0 + v) + g[2];
				y = g[3] * (u0 + u) + g[4] * (v0 + v) + g[5];
				x -= x0, y -= y0;
				if (x >= 0 && x < inWidth && y >= 0 && y < inHeight){
					Bicubic_hst(inPtr, channels*inWidth, channels, x, y,
						outPtr, channels*inHeight, u, v);
				}
				else{
					int i = INDEX(u, v, channels*outWidth, channels);
					outPtr[i] = outPtr[i + 1] = outPtr[i + 2] = 0;
				}
			}
		}
	}
}


void help()
{
	printf("\n\n----------- HUONG DAN -----------\n"
		"<filename>: Ten duong dan anh. D:\\img1.bmp\n"
		"[/f]: 6 tham so cua affice. [1 0 0 0 1 0]\n"
		"[/m]: phuong phap tai che mau. [NEARESTNEIGHBOR | BILINEAR | BICUBIC]\n"
		"[/s]: kich thuoc mau anh. [512 512]\n"
		"[/p]: lua chon che do. [CPU | GPU | OPENMP]\n"
		"[/t]: voi GPU, so threads 2d tren mot block. [8 8]\n"
		"[/b]: voi GPU, so blocks 2d grid. [16 16]\n");
}